home *** CD-ROM | disk | FTP | other *** search
/ AmigActive 21 / AACD 21.iso / AACD / Utilities / Ghostscript / src / gxpcopy.c < prev    next >
Encoding:
C/C++ Source or Header  |  2001-01-01  |  25.8 KB  |  896 lines
  1. /* Copyright (C) 1992, 2000 Aladdin Enterprises.  All rights reserved.
  2.   
  3.   This file is part of AFPL Ghostscript.
  4.   
  5.   AFPL Ghostscript is distributed with NO WARRANTY OF ANY KIND.  No author or
  6.   distributor accepts any responsibility for the consequences of using it, or
  7.   for whether it serves any particular purpose or works at all, unless he or
  8.   she says so in writing.  Refer to the Aladdin Free Public License (the
  9.   "License") for full details.
  10.   
  11.   Every copy of AFPL Ghostscript must include a copy of the License, normally
  12.   in a plain ASCII text file named PUBLIC.  The License grants you the right
  13.   to copy, modify and redistribute AFPL Ghostscript, but only under certain
  14.   conditions described in the License.  Among other things, the License
  15.   requires that the copyright notice and this notice be preserved on all
  16.   copies.
  17. */
  18.  
  19. /*$Id: gxpcopy.c,v 1.3 2000/09/19 19:00:39 lpd Exp $ */
  20. /* Path copying and flattening */
  21. #include "math_.h"
  22. #include "gx.h"
  23. #include "gserrors.h"
  24. #include "gconfigv.h"        /* for USE_FPU */
  25. #include "gxfixed.h"
  26. #include "gxfarith.h"
  27. #include "gxistate.h"        /* for access to line params */
  28. #include "gzpath.h"
  29.  
  30. /* Forward declarations */
  31. private void adjust_point_to_tangent(P3(segment *, const segment *,
  32.                     const gs_fixed_point *));
  33. private int monotonize_internal(P2(gx_path *, const curve_segment *));
  34.  
  35. /* Copy a path, optionally flattening or monotonizing it. */
  36. /* If the copy fails, free the new path. */
  37. int
  38. gx_path_copy_reducing(const gx_path *ppath_old, gx_path *ppath,
  39.               fixed fixed_flatness, const gs_imager_state *pis,
  40.               gx_path_copy_options options)
  41. {
  42.     const segment *pseg;
  43.     fixed flat = fixed_flatness;
  44.     gs_fixed_point expansion;
  45.     /*
  46.      * Since we're going to be adding to the path, unshare it
  47.      * before we start.
  48.      */
  49.     int code = gx_path_unshare(ppath);
  50.  
  51.     if (code < 0)
  52.     return code;
  53. #ifdef DEBUG
  54.     if (gs_debug_c('P'))
  55.     gx_dump_path(ppath_old, "before reducing");
  56. #endif
  57.     if (options & pco_for_stroke) {
  58.     /* Precompute the maximum expansion of the bounding box. */
  59.     double width = pis->line_params.half_width;
  60.  
  61.     expansion.x =
  62.         float2fixed((fabs(pis->ctm.xx) + fabs(pis->ctm.yx)) * width) * 2;
  63.     expansion.y =
  64.         float2fixed((fabs(pis->ctm.xy) + fabs(pis->ctm.yy)) * width) * 2;
  65.     }
  66.     pseg = (const segment *)(ppath_old->first_subpath);
  67.     while (pseg) {
  68.     switch (pseg->type) {
  69.         case s_start:
  70.         code = gx_path_add_point(ppath,
  71.                      pseg->pt.x, pseg->pt.y);
  72.         break;
  73.         case s_curve:
  74.         {
  75.             const curve_segment *pc = (const curve_segment *)pseg;
  76.  
  77.             if (fixed_flatness == max_fixed) {    /* don't flatten */
  78.             if (options & pco_monotonize)
  79.                 code = monotonize_internal(ppath, pc);
  80.             else
  81.                 code = gx_path_add_curve_notes(ppath,
  82.                      pc->p1.x, pc->p1.y, pc->p2.x, pc->p2.y,
  83.                        pc->pt.x, pc->pt.y, pseg->notes);
  84.             } else {
  85.             fixed x0 = ppath->position.x;
  86.             fixed y0 = ppath->position.y;
  87.             segment_notes notes = pseg->notes;
  88.             curve_segment cseg;
  89.             int k;
  90.  
  91.             if (options & pco_for_stroke) {
  92.                 /*
  93.                  * When flattening for stroking, the flatness
  94.                  * must apply to the outside of the resulting
  95.                  * stroked region.  We approximate this by
  96.                  * dividing the flatness by the ratio of the
  97.                  * expanded bounding box to the original
  98.                  * bounding box.  This is crude, but pretty
  99.                  * simple to calculate, and produces reasonably
  100.                  * good results.
  101.                  */
  102.                 fixed min01, max01, min23, max23;
  103.                 fixed ex, ey, flat_x, flat_y;
  104.  
  105. #define SET_EXTENT(r, c0, c1, c2, c3)\
  106.     BEGIN\
  107.     if (c0 < c1) min01 = c0, max01 = c1;\
  108.     else         min01 = c1, max01 = c0;\
  109.     if (c2 < c3) min23 = c2, max23 = c3;\
  110.     else         min23 = c3, max23 = c2;\
  111.     r = max(max01, max23) - min(min01, min23);\
  112.     END
  113.                 SET_EXTENT(ex, x0, pc->p1.x, pc->p2.x, pc->pt.x);
  114.                 SET_EXTENT(ey, y0, pc->p1.y, pc->p2.y, pc->pt.y);
  115. #undef SET_EXTENT
  116.                 /*
  117.                  * We check for the degenerate case specially
  118.                  * to avoid a division by zero.
  119.                  */
  120.                 if (ex == 0 || ey == 0)
  121.                 flat = 0;
  122.                 else {
  123.                 flat_x =
  124.                     fixed_mult_quo(fixed_flatness, ex,
  125.                            ex + expansion.x);
  126.                 flat_y =
  127.                     fixed_mult_quo(fixed_flatness, ey,
  128.                            ey + expansion.y);
  129.                 flat = min(flat_x, flat_y);
  130.                 }
  131.             }
  132.             k = gx_curve_log2_samples(x0, y0, pc, flat);
  133.             if (options & pco_accurate) {
  134.                 segment *start;
  135.                 segment *end;
  136.  
  137.                 /*
  138.                  * Add an extra line, which will become
  139.                  * the tangent segment.
  140.                  */
  141.                 code = gx_path_add_line_notes(ppath, x0, y0,
  142.                               notes);
  143.                 if (code < 0)
  144.                 break;
  145.                 start = ppath->current_subpath->last;
  146.                 notes |= sn_not_first;
  147.                 cseg = *pc;
  148.                 code = gx_flatten_sample(ppath, k, &cseg, notes);
  149.                 if (code < 0)
  150.                 break;
  151.                 /*
  152.                  * Adjust the first and last segments so that
  153.                  * they line up with the tangents.
  154.                  */
  155.                 end = ppath->current_subpath->last;
  156.                 if ((code = gx_path_add_line_notes(ppath,
  157.                               ppath->position.x,
  158.                               ppath->position.y,
  159.                         pseg->notes | sn_not_first)) < 0)
  160.                 break;
  161.                 adjust_point_to_tangent(start, start->next,
  162.                             &pc->p1);
  163.                 adjust_point_to_tangent(end, end->prev,
  164.                             &pc->p2);
  165.             } else {
  166.                 cseg = *pc;
  167.                 code = gx_flatten_sample(ppath, k, &cseg, notes);
  168.             }
  169.             }
  170.             break;
  171.         }
  172.         case s_line:
  173.         code = gx_path_add_line_notes(ppath,
  174.                        pseg->pt.x, pseg->pt.y, pseg->notes);
  175.         break;
  176.         case s_line_close:
  177.         code = gx_path_close_subpath(ppath);
  178.         break;
  179.         default:        /* can't happen */
  180.         code = gs_note_error(gs_error_unregistered);
  181.     }
  182.     if (code < 0) {
  183.         gx_path_new(ppath);
  184.         return code;
  185.     }
  186.     pseg = pseg->next;
  187.     }
  188.     if (path_last_is_moveto(ppath_old))
  189.     gx_path_add_point(ppath, ppath_old->position.x,
  190.               ppath_old->position.y);
  191.     if (ppath_old->bbox_set) {
  192.     if (ppath->bbox_set) {
  193.         ppath->bbox.p.x = min(ppath_old->bbox.p.x, ppath->bbox.p.x);
  194.         ppath->bbox.p.y = min(ppath_old->bbox.p.y, ppath->bbox.p.y);
  195.         ppath->bbox.q.x = max(ppath_old->bbox.q.x, ppath->bbox.q.x);
  196.         ppath->bbox.q.y = max(ppath_old->bbox.q.y, ppath->bbox.q.y);
  197.     } else {
  198.         ppath->bbox_set = true;
  199.         ppath->bbox = ppath_old->bbox;
  200.     }
  201.     }
  202. #ifdef DEBUG
  203.     if (gs_debug_c('P'))
  204.     gx_dump_path(ppath, "after reducing");
  205. #endif
  206.     return 0;
  207. }
  208.  
  209. /*
  210.  * Adjust one end of a line (the first or last line of a flattened curve)
  211.  * so it falls on the curve tangent.  The closest point on the line from
  212.  * (0,0) to (C,D) to a point (U,V) -- i.e., the point on the line at which
  213.  * a perpendicular line from the point intersects it -- is given by
  214.  *      T = (C*U + D*V) / (C^2 + D^2)
  215.  *      (X,Y) = (C*T,D*T)
  216.  * However, any smaller value of T will also work: the one we actually
  217.  * use is 0.25 * the value we just derived.  We must check that
  218.  * numerical instabilities don't lead to a negative value of T.
  219.  */
  220. private void
  221. adjust_point_to_tangent(segment * pseg, const segment * next,
  222.             const gs_fixed_point * p1)
  223. {
  224.     const fixed x0 = pseg->pt.x, y0 = pseg->pt.y;
  225.     const fixed fC = p1->x - x0, fD = p1->y - y0;
  226.  
  227.     /*
  228.      * By far the commonest case is that the end of the curve is
  229.      * horizontal or vertical.  Check for this specially, because
  230.      * we can handle it with far less work (and no floating point).
  231.      */
  232.     if (fC == 0) {
  233.     /* Vertical tangent. */
  234.     const fixed DT = arith_rshift(next->pt.y - y0, 2);
  235.  
  236.     if (fD == 0)
  237.         return;        /* anomalous case */
  238.     if_debug1('2', "[2]adjusting vertical: DT = %g\n",
  239.           fixed2float(DT));
  240.     if ((DT ^ fD) > 0)
  241.         pseg->pt.y = DT + y0;
  242.     } else if (fD == 0) {
  243.     /* Horizontal tangent. */
  244.     const fixed CT = arith_rshift(next->pt.x - x0, 2);
  245.  
  246.     if_debug1('2', "[2]adjusting horizontal: CT = %g\n",
  247.           fixed2float(CT));
  248.     if ((CT ^ fC) > 0)
  249.         pseg->pt.x = CT + x0;
  250.     } else {
  251.     /* General case. */
  252.     const double C = fC, D = fD;
  253.     const double T = (C * (next->pt.x - x0) + D * (next->pt.y - y0)) /
  254.     (C * C + D * D);
  255.  
  256.     if_debug3('2', "[2]adjusting: C = %g, D = %g, T = %g\n",
  257.           C, D, T);
  258.     if (T > 0) {
  259.         pseg->pt.x = arith_rshift((fixed) (C * T), 2) + x0;
  260.         pseg->pt.y = arith_rshift((fixed) (D * T), 2) + y0;
  261.     }
  262.     }
  263. }
  264.  
  265. /* ---------------- Curve flattening ---------------- */
  266.  
  267. #define x1 pc->p1.x
  268. #define y1 pc->p1.y
  269. #define x2 pc->p2.x
  270. #define y2 pc->p2.y
  271. #define x3 pc->pt.x
  272. #define y3 pc->pt.y
  273.  
  274. #ifdef DEBUG
  275. private void
  276. dprint_curve(const char *str, fixed x0, fixed y0, const curve_segment * pc)
  277. {
  278.     dlprintf9("%s p0=(%g,%g) p1=(%g,%g) p2=(%g,%g) p3=(%g,%g)\n",
  279.           str, fixed2float(x0), fixed2float(y0),
  280.           fixed2float(pc->p1.x), fixed2float(pc->p1.y),
  281.           fixed2float(pc->p2.x), fixed2float(pc->p2.y),
  282.           fixed2float(pc->pt.x), fixed2float(pc->pt.y));
  283. }
  284. #endif
  285.  
  286. /* Initialize a cursor for rasterizing a monotonic curve. */
  287. void
  288. gx_curve_cursor_init(curve_cursor * prc, fixed x0, fixed y0,
  289.              const curve_segment * pc, int k)
  290. {
  291.     fixed v01, v12;
  292.     int k2 = k + k, k3 = k2 + k;
  293.  
  294. #define bits_fit(v, n)\
  295.   (any_abs(v) <= max_fixed >> (n))
  296. /* The +2s are because of t3d and t2d, see below. */
  297. #define coeffs_fit(a, b, c)\
  298.   (k3 <= sizeof(fixed) * 8 - 3 &&\
  299.    bits_fit(a, k3 + 2) && bits_fit(b, k2 + 2) && bits_fit(c, k + 1))
  300.  
  301.     prc->k = k;
  302.     prc->p0.x = x0, prc->p0.y = y0;
  303.     prc->pc = pc;
  304.     /* Compute prc->a..c taking into account reversal of xy0/3 */
  305.     /* in curve_x_at_y. */
  306.     {
  307.     fixed w0, w1, w2, w3;
  308.  
  309.     if (y0 < y3)
  310.         w0 = x0, w1 = x1, w2 = x2, w3 = x3;
  311.     else
  312.         w0 = x3, w1 = x2, w2 = x1, w3 = x0;
  313.     curve_points_to_coefficients(w0, w1, w2, w3,
  314.                      prc->a, prc->b, prc->c, v01, v12);
  315.     }
  316.     prc->double_set = false;
  317.     prc->fixed_limit =
  318.     (coeffs_fit(prc->a, prc->b, prc->c) ? (1 << k) - 1 : -1);
  319.     /* Initialize the cache. */
  320.     prc->cache.ky0 = prc->cache.ky3 = y0;
  321.     prc->cache.xl = x0;
  322.     prc->cache.xd = 0;
  323. }
  324.  
  325. /*
  326.  * Determine the X value on a monotonic curve at a given Y value.  It turns
  327.  * out that we use so few points on the curve that it's actually faster to
  328.  * locate the desired point by recursive subdivision each time than to try
  329.  * to keep a cursor that we move incrementally.  What's even more surprising
  330.  * is that if floating point arithmetic is reasonably fast, it's faster to
  331.  * compute the X value at the desired point explicitly than to do the
  332.  * recursive subdivision on X as well as Y.
  333.  */
  334. #define SUBDIVIDE_X USE_FPU_FIXED
  335. fixed
  336. gx_curve_x_at_y(curve_cursor * prc, fixed y)
  337. {
  338.     fixed xl, xd;
  339.     fixed yd, yrel;
  340.  
  341.     /* Check the cache before doing anything else. */
  342.     if (y >= prc->cache.ky0 && y <= prc->cache.ky3) {
  343.     yd = prc->cache.ky3 - prc->cache.ky0;
  344.     yrel = y - prc->cache.ky0;
  345.     xl = prc->cache.xl;
  346.     xd = prc->cache.xd;
  347.     goto done;
  348.     } {
  349. #define x0 prc->p0.x
  350. #define y0 prc->p0.y
  351.     const curve_segment *pc = prc->pc;
  352.  
  353.     /* Reduce case testing by ensuring y3 >= y0. */
  354.     fixed cy0 = y0, cy1, cy2, cy3 = y3;
  355.     fixed cx0;
  356.  
  357. #if SUBDIVIDE_X
  358.     fixed cx1, cx2, cx3;
  359.  
  360. #else
  361.     int t = 0;
  362.  
  363. #endif
  364.     int k, i;
  365.  
  366.     if (cy0 > cy3)
  367.         cx0 = x3,
  368. #if SUBDIVIDE_X
  369.         cx1 = x2, cx2 = x1, cx3 = x0,
  370. #endif
  371.         cy0 = y3, cy1 = y2, cy2 = y1, cy3 = y0;
  372.     else
  373.         cx0 = x0,
  374. #if SUBDIVIDE_X
  375.         cx1 = x1, cx2 = x2, cx3 = x3,
  376. #endif
  377.         cy1 = y1, cy2 = y2;
  378. #define midpoint_fast(a,b)\
  379.   arith_rshift_1((a) + (b) + 1)
  380.     for (i = k = prc->k; i > 0; --i) {
  381.         fixed ym = midpoint_fast(cy1, cy2);
  382.         fixed yn = ym + arith_rshift(cy0 - cy1 - cy2 + cy3 + 4, 3);
  383.  
  384. #if SUBDIVIDE_X
  385.         fixed xm = midpoint_fast(cx1, cx2);
  386.         fixed xn = xm + arith_rshift(cx0 - cx1 - cx2 + cx3 + 4, 3);
  387.  
  388. #else
  389.         t <<= 1;
  390. #endif
  391.  
  392.         if (y < yn)
  393. #if SUBDIVIDE_X
  394.         cx1 = midpoint_fast(cx0, cx1),
  395.             cx2 = midpoint_fast(cx1, xm),
  396.             cx3 = xn,
  397. #endif
  398.             cy1 = midpoint_fast(cy0, cy1),
  399.             cy2 = midpoint_fast(cy1, ym),
  400.             cy3 = yn;
  401.         else
  402. #if SUBDIVIDE_X
  403.         cx2 = midpoint_fast(cx2, cx3),
  404.             cx1 = midpoint_fast(xm, cx2),
  405.             cx0 = xn,
  406. #else
  407.         t++,
  408. #endif
  409.             cy2 = midpoint_fast(cy2, cy3),
  410.             cy1 = midpoint_fast(ym, cy2),
  411.             cy0 = yn;
  412.     }
  413. #if SUBDIVIDE_X
  414.     xl = cx0;
  415.     xd = cx3 - cx0;
  416. #else
  417.     {
  418.         fixed a = prc->a, b = prc->b, c = prc->c;
  419.  
  420. /* We must use (1 << k) >> 1 instead of 1 << (k - 1) in case k == 0. */
  421. #define compute_fixed(a, b, c)\
  422.   arith_rshift(arith_rshift(arith_rshift(a * t3, k) + b * t2, k)\
  423.            + c * t + ((1 << k) >> 1), k)
  424. #define compute_diff_fixed(a, b, c)\
  425.   arith_rshift(arith_rshift(arith_rshift(a * t3d, k) + b * t2d, k)\
  426.            + c, k)
  427.  
  428.         /* use multiply if possible */
  429. #define np2(n) (1.0 / (1L << (n)))
  430.         static const double k_denom[11] =
  431.         {np2(0), np2(1), np2(2), np2(3), np2(4),
  432.          np2(5), np2(6), np2(7), np2(8), np2(9), np2(10)
  433.         };
  434.         static const double k2_denom[11] =
  435.         {np2(0), np2(2), np2(4), np2(6), np2(8),
  436.          np2(10), np2(12), np2(14), np2(16), np2(18), np2(20)
  437.         };
  438.         static const double k3_denom[11] =
  439.         {np2(0), np2(3), np2(6), np2(9), np2(12),
  440.          np2(15), np2(18), np2(21), np2(24), np2(27), np2(30)
  441.         };
  442.         double den1, den2;
  443.  
  444. #undef np2
  445.  
  446. #define setup_floating(da, db, dc, a, b, c)\
  447.   (k >= countof(k_denom) ?\
  448.    (den1 = ldexp(1.0, -k),\
  449.     den2 = den1 * den1,\
  450.     da = (den2 * den1) * a,\
  451.     db = den2 * b,\
  452.     dc = den1 * c) :\
  453.    (da = k3_denom[k] * a,\
  454.     db = k2_denom[k] * b,\
  455.     dc = k_denom[k] * c))
  456. #define compute_floating(da, db, dc)\
  457.   ((fixed)(da * t3 + db * t2 + dc * t + 0.5))
  458. #define compute_diff_floating(da, db, dc)\
  459.   ((fixed)(da * t3d + db * t2d + dc))
  460.  
  461.         if (t <= prc->fixed_limit) {    /* We can do everything in integer/fixed point. */
  462.         int t2 = t * t, t3 = t2 * t;
  463.         int t3d = (t2 + t) * 3 + 1, t2d = t + t + 1;
  464.  
  465.         xl = compute_fixed(a, b, c) + cx0;
  466.         xd = compute_diff_fixed(a, b, c);
  467. #ifdef DEBUG
  468.         {
  469.             double fa, fb, fc;
  470.             fixed xlf, xdf;
  471.  
  472.             setup_floating(fa, fb, fc, a, b, c);
  473.             xlf = compute_floating(fa, fb, fc) + cx0;
  474.             xdf = compute_diff_floating(fa, fb, fc);
  475.             if (any_abs(xlf - xl) > fixed_epsilon ||
  476.             any_abs(xdf - xd) > fixed_epsilon
  477.             )
  478.             dlprintf9("Curve points differ: k=%d t=%d a,b,c=%g,%g,%g\n   xl,xd fixed=%g,%g floating=%g,%g\n",
  479.                   k, t,
  480.                  fixed2float(a), fixed2float(b), fixed2float(c),
  481.                   fixed2float(xl), fixed2float(xd),
  482.                   fixed2float(xlf), fixed2float(xdf));
  483. /*xl = xlf, xd = xdf; */
  484.         }
  485. #endif
  486.         } else {        /*
  487.                  * Either t3 (and maybe t2) won't fit in an int, or more
  488.                  * likely the result of the multiplies won't fit.
  489.                  */
  490. #define fa prc->da
  491. #define fb prc->db
  492. #define fc prc->dc
  493.         if (!prc->double_set) {
  494.             setup_floating(fa, fb, fc, a, b, c);
  495.             prc->double_set = true;
  496.         }
  497.         if (t < 1L << ((sizeof(long) * 8 - 1) / 3)) {    /*
  498.                                  * t3 (and maybe t2) might not fit in an int, but they
  499.                                  * will fit in a long.  If we have slow floating point,
  500.                                  * do the computation in double-precision fixed point,
  501.                                  * otherwise do it in fixed point.
  502.                                  */
  503.             long t2 = (long)t * t, t3 = t2 * t;
  504.             long t3d = (t2 + t) * 3 + 1, t2d = t + t + 1;
  505.  
  506.             xl = compute_floating(fa, fb, fc) + cx0;
  507.             xd = compute_diff_floating(fa, fb, fc);
  508.         } else {    /*
  509.                  * t3 (and maybe t2) don't even fit in a long.
  510.                  * Do the entire computation in floating point.
  511.                  */
  512.             double t2 = (double)t * t, t3 = t2 * t;
  513.             double t3d = (t2 + t) * 3 + 1, t2d = t + t + 1;
  514.  
  515.             xl = compute_floating(fa, fb, fc) + cx0;
  516.             xd = compute_diff_floating(fa, fb, fc);
  517.         }
  518. #undef fa
  519. #undef fb
  520. #undef fc
  521.         }
  522.     }
  523. #endif /* (!)SUBDIVIDE_X */
  524.  
  525.     /* Update the cache. */
  526.     prc->cache.ky0 = cy0;
  527.     prc->cache.ky3 = cy3;
  528.     prc->cache.xl = xl;
  529.     prc->cache.xd = xd;
  530.     yd = cy3 - cy0;
  531.     yrel = y - cy0;
  532. #undef x0
  533. #undef y0
  534.     }
  535.   done:
  536.     /*
  537.      * Now interpolate linearly between current and next.
  538.      * We know that 0 <= yrel < yd.
  539.      * It's unlikely but possible that cy0 = y = cy3:
  540.      * handle this case specially.
  541.      */
  542.     if (yrel == 0)
  543.     return xl;
  544.     /*
  545.      * Compute in fixed point if possible.
  546.      */
  547. #define HALF_FIXED_BITS ((fixed)1 << (sizeof(fixed) * 4))
  548.     if (yrel < HALF_FIXED_BITS) {
  549.     if (xd >= 0) {
  550.         if (xd < HALF_FIXED_BITS)
  551.         return (ufixed)xd * (ufixed)yrel / (ufixed)yd + xl;
  552.     } else {
  553.         if (xd > -HALF_FIXED_BITS) {
  554.         /* Be careful to take the floor of the result. */
  555.         ufixed num = (ufixed)(-xd) * (ufixed)yrel;
  556.         ufixed quo = num / (ufixed)yd;
  557.  
  558.         if (quo * (ufixed)yd != num)
  559.             quo += fixed_epsilon;
  560.         return xl - (fixed)quo;
  561.         }
  562.     }
  563.     }
  564. #undef HALF_FIXED_BITS
  565.     return fixed_mult_quo(xd, yrel, yd) + xl;
  566. }
  567.  
  568. #undef x1
  569. #undef y1
  570. #undef x2
  571. #undef y2
  572. #undef x3
  573. #undef y3
  574.  
  575. /* ---------------- Monotonic curves ---------------- */
  576.  
  577. /* Test whether a path is free of non-monotonic curves. */
  578. bool
  579. gx_path_is_monotonic(const gx_path * ppath)
  580. {
  581.     const segment *pseg = (const segment *)(ppath->first_subpath);
  582.     gs_fixed_point pt0;
  583.  
  584.     while (pseg) {
  585.     switch (pseg->type) {
  586.         case s_start:
  587.         {
  588.             const subpath *psub = (const subpath *)pseg;
  589.  
  590.             /* Skip subpaths without curves. */
  591.             if (!psub->curve_count)
  592.             pseg = psub->last;
  593.         }
  594.         break;
  595.         case s_curve:
  596.         {
  597.             const curve_segment *pc = (const curve_segment *)pseg;
  598.             double t[2];
  599.             int nz = gx_curve_monotonic_points(pt0.y,
  600.                        pc->p1.y, pc->p2.y, pc->pt.y, t);
  601.  
  602.             if (nz != 0)
  603.             return false;
  604.             nz = gx_curve_monotonic_points(pt0.x,
  605.                        pc->p1.x, pc->p2.x, pc->pt.x, t);
  606.             if (nz != 0)
  607.             return false;
  608.         }
  609.         break;
  610.         default:
  611.         ;
  612.     }
  613.     pt0 = pseg->pt;
  614.     pseg = pseg->next;
  615.     }
  616.     return true;
  617. }
  618.  
  619. /* Monotonize a curve, by splitting it if necessary. */
  620. /* In the worst case, this could split the curve into 9 pieces. */
  621. private int
  622. monotonize_internal(gx_path * ppath, const curve_segment * pc)
  623. {
  624.     fixed x0 = ppath->position.x, y0 = ppath->position.y;
  625.     segment_notes notes = pc->notes;
  626.     double t[2];
  627.  
  628. #define max_segs 9
  629.     curve_segment cs[max_segs];
  630.     const curve_segment *pcs;
  631.     curve_segment *pcd;
  632.     int i, j, nseg;
  633.     int nz;
  634.  
  635.     /* Monotonize in Y. */
  636.     nz = gx_curve_monotonic_points(y0, pc->p1.y, pc->p2.y, pc->pt.y, t);
  637.     nseg = max_segs - 1 - nz;
  638.     pcd = cs + nseg;
  639.     if (nz == 0)
  640.     *pcd = *pc;
  641.     else {
  642.     gx_curve_split(x0, y0, pc, t[0], pcd, pcd + 1);
  643.     if (nz == 2)
  644.         gx_curve_split(pcd->pt.x, pcd->pt.y, pcd + 1,
  645.                (t[1] - t[0]) / (1 - t[0]),
  646.                pcd + 1, pcd + 2);
  647.     }
  648.  
  649.     /* Monotonize in X. */
  650.     for (pcs = pcd, pcd = cs, j = nseg; j < max_segs; ++pcs, ++j) {
  651.     nz = gx_curve_monotonic_points(x0, pcs->p1.x, pcs->p2.x,
  652.                        pcs->pt.x, t);
  653.  
  654.     if (nz == 0)
  655.         *pcd = *pcs;
  656.     else {
  657.         gx_curve_split(x0, y0, pcs, t[0], pcd, pcd + 1);
  658.         if (nz == 2)
  659.         gx_curve_split(pcd->pt.x, pcd->pt.y, pcd + 1,
  660.                    (t[1] - t[0]) / (1 - t[0]),
  661.                    pcd + 1, pcd + 2);
  662.     }
  663.     pcd += nz + 1;
  664.     x0 = pcd[-1].pt.x;
  665.     y0 = pcd[-1].pt.y;
  666.     }
  667.     nseg = pcd - cs;
  668.  
  669.     /* Add the segment(s) to the output. */
  670. #ifdef DEBUG
  671.     if (gs_debug_c('2')) {
  672.     int pi;
  673.     gs_fixed_point pp0;
  674.  
  675.     pp0 = ppath->position;
  676.     if (nseg == 1)
  677.         dprint_curve("[2]No split", pp0.x, pp0.y, pc);
  678.     else {
  679.         dlprintf1("[2]Split into %d segments:\n", nseg);
  680.         dprint_curve("[2]Original", pp0.x, pp0.y, pc);
  681.         for (pi = 0; pi < nseg; ++pi) {
  682.         dprint_curve("[2] =>", pp0.x, pp0.y, cs + pi);
  683.         pp0 = cs[pi].pt;
  684.         }
  685.     }
  686.     }
  687. #endif
  688.     for (pcs = cs, i = 0; i < nseg; ++pcs, ++i) {
  689.     int code = gx_path_add_curve_notes(ppath, pcs->p1.x, pcs->p1.y,
  690.                        pcs->p2.x, pcs->p2.y,
  691.                        pcs->pt.x, pcs->pt.y,
  692.                        notes |
  693.                        (i > 0 ? sn_not_first :
  694.                         sn_none));
  695.  
  696.     if (code < 0)
  697.         return code;
  698.     }
  699.  
  700.     return 0;
  701. }
  702.  
  703. /*
  704.  * Split a curve if necessary into pieces that are monotonic in X or Y as a
  705.  * function of the curve parameter t.  This allows us to rasterize curves
  706.  * directly without pre-flattening.  This takes a fair amount of analysis....
  707.  * Store the values of t of the split points in pst[0] and pst[1].  Return
  708.  * the number of split points (0, 1, or 2).
  709.  */
  710. int
  711. gx_curve_monotonic_points(fixed v0, fixed v1, fixed v2, fixed v3,
  712.               double pst[2])
  713. {
  714.     /*
  715.        Let
  716.        v(t) = a*t^3 + b*t^2 + c*t + d, 0 <= t <= 1.
  717.        Then
  718.        dv(t) = 3*a*t^2 + 2*b*t + c.
  719.        v(t) has a local minimum or maximum (or inflection point)
  720.        precisely where dv(t) = 0.  Now the roots of dv(t) = 0 (i.e.,
  721.        the zeros of dv(t)) are at
  722.        t =  ( -2*b +/- sqrt(4*b^2 - 12*a*c) ) / 6*a
  723.        =    ( -b +/- sqrt(b^2 - 3*a*c) ) / 3*a
  724.        (Note that real roots exist iff b^2 >= 3*a*c.)
  725.        We want to know if these lie in the range (0..1).
  726.        (The endpoints don't count.)  Call such a root a "valid zero."
  727.        Since computing the roots is expensive, we would like to have
  728.        some cheap tests to filter out cases where they don't exist
  729.        (i.e., where the curve is already monotonic).
  730.      */
  731.     fixed v01, v12, a, b, c, b2, a3;
  732.     fixed dv_end, b2abs, a3abs;
  733.  
  734.     curve_points_to_coefficients(v0, v1, v2, v3, a, b, c, v01, v12);
  735.     b2 = b << 1;
  736.     a3 = (a << 1) + a;
  737.     /*
  738.        If a = 0, the only possible zero is t = -c / 2*b.
  739.        This zero is valid iff sign(c) != sign(b) and 0 < |c| < 2*|b|.
  740.      */
  741.     if (a == 0) {
  742.     if ((b ^ c) < 0 && any_abs(c) < any_abs(b2) && c != 0) {
  743.         *pst = (double)(-c) / b2;
  744.         return 1;
  745.     } else
  746.         return 0;
  747.     }
  748.     /*
  749.        Iff a curve is horizontal at t = 0, c = 0.  In this case,
  750.        there can be at most one other zero, at -2*b / 3*a.
  751.        This zero is valid iff sign(a) != sign(b) and 0 < 2*|b| < 3*|a|.
  752.      */
  753.     if (c == 0) {
  754.     if ((a ^ b) < 0 && any_abs(b2) < any_abs(a3) && b != 0) {
  755.         *pst = (double)(-b2) / a3;
  756.         return 1;
  757.     } else
  758.         return 0;
  759.     }
  760.     /*
  761.        Similarly, iff a curve is horizontal at t = 1, 3*a + 2*b + c = 0.
  762.        In this case, there can be at most one other zero,
  763.        at -1 - 2*b / 3*a, iff sign(a) != sign(b) and 1 < -2*b / 3*a < 2,
  764.        i.e., 3*|a| < 2*|b| < 6*|a|.
  765.      */
  766.     else if ((dv_end = a3 + b2 + c) == 0) {
  767.     if ((a ^ b) < 0 &&
  768.         (b2abs = any_abs(b2)) > (a3abs = any_abs(a3)) &&
  769.         b2abs < a3abs << 1
  770.         ) {
  771.         *pst = (double)(-b2 - a3) / a3;
  772.         return 1;
  773.     } else
  774.         return 0;
  775.     }
  776.     /*
  777.        If sign(dv_end) != sign(c), at least one valid zero exists,
  778.        since dv(0) and dv(1) have opposite signs and hence
  779.        dv(t) must be zero somewhere in the interval [0..1].
  780.      */
  781.     else if ((dv_end ^ c) < 0);
  782.     /*
  783.        If sign(a) = sign(b), no valid zero exists,
  784.        since dv is monotonic on [0..1] and has the same sign
  785.        at both endpoints.
  786.      */
  787.     else if ((a ^ b) >= 0)
  788.     return 0;
  789.     /*
  790.        Otherwise, dv(t) may be non-monotonic on [0..1]; it has valid zeros
  791.        iff its sign anywhere in this interval is different from its sign
  792.        at the endpoints, which occurs iff it has an extremum in this
  793.        interval and the extremum is of the opposite sign from c.
  794.        To find this out, we look for the local extremum of dv(t)
  795.        by observing
  796.        d2v(t) = 6*a*t + 2*b
  797.        which has a zero only at
  798.        t1 = -b / 3*a
  799.        Now if t1 <= 0 or t1 >= 1, no valid zero exists.
  800.        Note that we just determined that sign(a) != sign(b), so we know t1 > 0.
  801.      */
  802.     else if (any_abs(b) >= any_abs(a3))
  803.     return 0;
  804.     /*
  805.        Otherwise, we just go ahead with the computation of the roots,
  806.        and test them for being in the correct range.  Note that a valid
  807.        zero is an inflection point of v(t) iff d2v(t) = 0; we don't
  808.        bother to check for this case, since it's rare.
  809.      */
  810.     {
  811.     double nbf = (double)(-b);
  812.     double a3f = (double)a3;
  813.     double radicand = nbf * nbf - a3f * c;
  814.  
  815.     if (radicand < 0) {
  816.         if_debug1('2', "[2]negative radicand = %g\n", radicand);
  817.         return 0;
  818.     } {
  819.         double root = sqrt(radicand);
  820.         int nzeros = 0;
  821.         double z = (nbf - root) / a3f;
  822.  
  823.         /*
  824.          * We need to return the zeros in the correct order.
  825.          * We know that root is non-negative, but a3f may be either
  826.          * positive or negative, so we need to check the ordering
  827.          * explicitly.
  828.          */
  829.         if_debug2('2', "[2]zeros at %g, %g\n", z, (nbf + root) / a3f);
  830.         if (z > 0 && z < 1)
  831.         *pst = z, nzeros = 1;
  832.         if (root != 0) {
  833.         z = (nbf + root) / a3f;
  834.         if (z > 0 && z < 1) {
  835.             if (nzeros && a3f < 0)    /* order is reversed */
  836.             pst[1] = *pst, *pst = z;
  837.             else
  838.             pst[nzeros] = z;
  839.             nzeros++;
  840.         }
  841.         }
  842.         return nzeros;
  843.     }
  844.     }
  845. }
  846.  
  847. /*
  848.  * Split a curve at an arbitrary point t.  The above midpoint split is a
  849.  * special case of this with t = 0.5.
  850.  */
  851. void
  852. gx_curve_split(fixed x0, fixed y0, const curve_segment * pc, double t,
  853.            curve_segment * pc1, curve_segment * pc2)
  854. {                /*
  855.                  * If the original function was v(t), we want to compute the points
  856.                  * for the functions v1(T) = v(t * T) and v2(T) = v(t + (1 - t) * T).
  857.                  * Straightforwardly,
  858.                  *      v1(T) = a*t^3*T^3 + b*t^2*T^2 + c*t*T + d
  859.                  * i.e.
  860.                  *      a1 = a*t^3, b1 = b*t^2, c1 = c*t, d1 = d.
  861.                  * Similarly,
  862.                  *      v2(T) = a*[t + (1-t)*T]^3 + b*[t + (1-t)*T]^2 +
  863.                  *              c*[t + (1-t)*T] + d
  864.                  *            = a*[(1-t)^3*T^3 + 3*t*(1-t)^2*T^2 + 3*t^2*(1-t)*T +
  865.                  *                 t^3] + b*[(1-t)^2*T^2 + 2*t*(1-t)*T + t^2] +
  866.                  *                 c*[(1-t)*T + t] + d
  867.                  *            = a*(1-t)^3*T^3 + [a*3*t + b]*(1-t)^2*T^2 +
  868.                  *                 [a*3*t^2 + b*2*t + c]*(1-t)*T +
  869.                  *                 a*t^3 + b*t^2 + c*t + d
  870.                  * We do this in the simplest way, namely, we convert the points to
  871.                  * coefficients, do the arithmetic, and convert back.  It would
  872.                  * obviously be faster to do the arithmetic directly on the points,
  873.                  * as the midpoint code does; this is just an implementation issue
  874.                  * that we can revisit if necessary.
  875.                  */
  876.     double t2 = t * t, t3 = t2 * t;
  877.     double omt = 1 - t, omt2 = omt * omt, omt3 = omt2 * omt;
  878.     fixed v01, v12, a, b, c, na, nb, nc;
  879.  
  880.     if_debug1('2', "[2]splitting at t = %g\n", t);
  881. #define compute_seg(v0, v)\
  882.     curve_points_to_coefficients(v0, pc->p1.v, pc->p2.v, pc->pt.v,\
  883.                      a, b, c, v01, v12);\
  884.     na = (fixed)(a * t3), nb = (fixed)(b * t2), nc = (fixed)(c * t);\
  885.     curve_coefficients_to_points(na, nb, nc, v0,\
  886.                      pc1->p1.v, pc1->p2.v, pc1->pt.v);\
  887.     na = (fixed)(a * omt3);\
  888.     nb = (fixed)((a * t * 3 + b) * omt2);\
  889.     nc = (fixed)((a * t2 * 3 + b * 2 * t + c) * omt);\
  890.     curve_coefficients_to_points(na, nb, nc, pc1->pt.v,\
  891.                      pc2->p1.v, pc2->p2.v, pc2->pt.v)
  892.     compute_seg(x0, x);
  893.     compute_seg(y0, y);
  894. #undef compute_seg
  895. }
  896.